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Abstract 

The application of concepts from equilibrium statistical mechanics to out-of-equilibrium systems 
has a long history of describing diverse systems ranging from glasses to granular materials. For 
dissipative jammed systems- particulate grains or droplets- a key concept is to replace the en- 
ergy ensemble describing conservative systems by the volume-stress ensemble. Here, we test the 
applicability of the volume-stress ensemble to describe the jamming transition by comparing the 
jammed configurations obtained by dynamics with those averaged over the ensemble as a probe 
of ergodicity. Agreement between both methods suggests the idea of "thermalization" at a given 
angoricity and compactivity. We elucidate the thermodynamic order of the jamming transition by 
showing the absence of critical fluctuations in static observables like pressure and volume. The 
approach allows to calculate observables such as the entropy, volume, pressure, coordination num- 
ber and distribution of forces to characterize the scaling laws near the jamming transition from a 
statistical mechanics viewpoint. 
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A granular system compresses into a mechanically stable configuration at a nonzero pres- 



sure in response to the application of an external strain 



-|3[. This process is typically 



referred to as the jamming transition and occurs at a critical volume fraction 0^ jsl. The ap- 
plication of a subsequent external pressure with the concomitant particle rearrangements and 
compression results in a set of configurations characterized by the system volume V = NVg / 
(0 is the volume fraction of particles of volume Vg) and applied external stress or pressure 
p (for simplicity we assume isotropic states). It has been long argued whether the jamming 
transition is a first-order transition at the discontinuity in the average coordination number, 
(Z), or a second-order transition with the power- law scaling of the system's pressure as the 
system approaches jamming with — 0^ ~^ O"*" jj-?]. Previous work js-lOj has proposed to 
explain the jamming transition by a field theory in the pressure ensemble. Here, we use the 
idea of "thermalization" of an ensemble of mechanically stable granular materials at a given 
volume and pressure to study the jamming transition from a thermodynamic viewpoint. 

For a fixed number of grains, there exist many jammed states Uj confined by the ex- 
ternal pressure p in a volume V. In an effort to describe the nature of this nonequilibrium 
system from a statistical mechanics perspective, a pressure-volume ensemble jsl. [l2ll5| was 
introduced for jammed matter. In the canonical the probability of a state is given by 
exp[—W{dS/dV) — T{dS/dT)], where S is the entropy of the system, W is the volume 
function measuring the volume of the system as a function of the particle coordinates and 
r = pV is the boundary stress (or internal virial) js] of the system. Just as dE/dS = T is 
the temperature in equilibrium system, the temperature-like variables in jammed systems 
are the compactivity X = dV/dS [12| and the angoricity A = dT/dS |l3 |. 

In a recent paper 16(| the compactivity was used to describe frictional hard spheres in 
the volume ensemble. Here, we test the validity of the statistical approach in the combined 
pressure-volume ensemble to describe deformable, frictionless particles, such as emulsion 
systems jammed under osmotic pressure near the jamming transition [17|. We demonstrate 
that the jamming transition can be probed thermodynamically by the angoricity A and 
the compactivity X. The calculation of jamming "heat" capacities characterizes the sys- 
tem fluctuations and shows the lack of critical fluctuations in the static quantities as the 
jamming transition point is approached from above — )■ 0+. Thus, the thermodynami- 
cal viewpoint determines the order of the phase transition and allows one to calculate the 
physical observables near jamming. 
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I. RESULTS 



In general, if the density of states g(T,(f)) in the space of jammed configurations (defined 
as the probabihty of finding a jammed state at a given (F, (p) ai A = oo) is known, then 
calculations of macroscopic observables, like pressure p and average coordination number Z 



as a function of (J), can be performed by the canonical ensemble average [o]. 



1^ at a given 



volume 
and 



{p{a, 0))ens = 1 1^ P 0) e-"^ dr, (1) 



1 

where the canonical partition function is Z = g(T, 0)e~"'"dr and the density of states is 



normalized as g(T,(j))dT = 1. The inverse angoricity is defined as a = 1/A = dS/dT. 

At the jamming transition the system reaches isostatic equilibrium, such that the stresses 
are exactly balanced in the resulting configuration, and there exists a unique solution to 
the interparticle force equations satisfjang mechanical equilibrium. It is well known that 
observables present power-law scaling 



ing met 

m 



(P)dyn ~ (0 - (f>cr , (^)dyn - Z, ^ {(j) - (3) 

where a = 3/2 and b = 1 /2 for Hertzian spheres and = 6 is the coordination number at the 
isostatic point (J-point) [4]. The average (■ ■ ■ )dyn indicates that these quantities are obtained 
by averaging over packings generated dynamically in either simulations or experiments as 
opposed to the ensemble average over configurations (■ ■ ■ )ens of Eqs. ([I])-([2]). Comparing 
the ensemble calculations, Eq. ([I])-(l2]), with the direct dynamical measurements, Eq. (jSI), 
provides a basic test of the ergodic hypothesis for the statistical ensemble. 

Our approach is the following: We first perform an exhaustive enumeration of configura- 
tions to calculate g(T,(f)) and obtain (p(a,0))ens as a function of a for a given using Eq. 
([1]). Then, we obtain the angoricity by comparing the pressure in the ensemble average with 
the one obtained following the dynamical evolution with Molecular Dynamics (MD) simula- 
tions. By setting {p{a, 0))cns = (p)dyn, we obtain the angoricity as a function of 0. By virtue 
of obtaining a{(f)), all the other observables can be calculated in the ensemble formulation. 
The ultimate test of ergodicity is realized by comparing the remaining ensemble observables 
with the corresponding direct dynamical measures. 



Ensemble calculations. — The density of jammed states g{r,(f)) is calculated in the 
ramework of the potential energy landscape (PEL) formulation introduced by Goldstein 



18| and Stillinger- Weber [19|, |20| to describe supercooled liquids. In the case of frictionless 



jammed systems, the mechanically stable configurations are defined as the local minima of 
the potential energy surface (PES) of the system {4, 11| (see Fig. [T] inset for a schematic 
representation). In the simulations, two spherical soft particles in contact interact via a 
normal Hertz force m,F„. (,.)^ where . the .o^a. overlap between the spheres 
under deformation and 5 = 1.5, in a periodically repeated cube [the interparticle potential 
energy is oc {5r)^^^, see Materials and Methods Section [III A] . The Hertz potential is 
chosen for its general applicability to granular materials. The results are expected to be 



independent of the form of the potential. Details of the algorithms |22l . |23| to find the local 
minima of the PES (zero-order saddles) are in the Materials and Methods Section IIII B[ 
Figure [1] shows g{r,(j)) versus F for different volume fractions. 

MD calculations. — The pressure (p)dyn as a function of is calculated by performing 
MD simulations. Packings are prepared by compressing a gas of particles from an initial 
(unjammed) low volume fraction to a final jammed state. This procedure simulates a dy- 



24j : details appear in Materials and Methods Section [HI CI 



namical packing preparation 
We obtain (Fig. [2^) 

(p)dyn=P0 {<!> - <f>cy-'' , (4) 

where 0c = 0.6077 is the volume fraction corresponding to the isostatic point J j4| following 
Eq. ([3]) and po = lO.SMPa. This critical value and the exponent, a = 1.65, are slightly 
different than the values obtained for larger systems (a = (5) [4] • However, our purpose is to 
use the same system in the dynamical calculation and the exact enumeration for a proper 
comparison. 

Calculation of angoricity. — For each we use g{r,(f)) to calculate (p(a))ens by Eq. 
dl]). Then, we obtain a(0) by setting (p(a;,0))ens = (p)dyn for every (see Figs. [^ and ITU] 
and Materials and Methods Section [IIID|) . The resulting equation of state a(0) is plotted 
in Fig. |2j3 and shows that the angoricity follows a power-law, near 0^, of the form: 



A oc 



(5) 



where the angoricity exponent is 7 = 2.5. The result is consistent with 7 = 5 + 1.0, 
suggesting that A (x V (x F^r. Angoricity is a measure of the number of ways the stress 



can be distributed in a given volume. Since the stresses have a unique solution for a given 
configuration at the isostatic point, (pc-, the corresponding angoricity vanishes. At higher 
pressure, the system is determined by multiple degrees of freedom satisfying mechanical 
equilibrium, leading to a higher stress temperature, A. The angoricity can also be viewed as 
a scale of stability for the system at different volume fractions. Systems jammed at larger 
volume fractions require higher angoricity (higher driving force) to rearrange. 

Test of ergodicity. — In principle, using the inverse angoricity, a, from Eq. ([5]) we can 
calculate any macroscopic statistical observable (-B)ens at a given volume by performing the 
ensemble average jlol |: 

We test the ergodic hypothesis in the Edwards's ensemble by comparing Eq. ([H]) with the 
corresponding value obtained with MD simulations averaged over (250) sample packings, 
generated dynamically: 

250 

(5(0))dyn=— 5^5. (7) 
i=l 

The comparison is realized by measuring the average coordination number, (Z), the 
average force and the distribution of interparticle forces. We calculate {Z)cns by Eq. ([2]) and 
(Z)dyn as in Eq. ([7]). Figures [3JA. and[3]B show that the two independent estimations of the 
coordination number agree very well: {Z)ens = {Z)dyn- 

We calculate the ensemble average force (-F)ons and the average over all the MD packings, 
{F)dyn and find that they coincide very closely (see Fig. |3p). The full distribution of 
inter-particle forces for jammed systems is also an important observable which has been 



extensively studied in previous works 



25 



261]. The force distribution is calculated in the 



ensemble Pens{F/F) by averaging the force distribution for every configuration in the PES 
(see Materials and Methods Section UlI D|) . Figure [3p shows the distribution functions. The 



peak of the distribution shown in Fig. [3p indicates that the systems are jammed j^, |25|, |26 |. 
Besides the exact shape of the distribution, the similarity between the ensemble and the 
dynamical calculations shown in Fig. |3p is significant. The study of (Z), (F) and P{F/F) 
reveals that the statistical ensemble can predict the macroscopic observables obtained in 
MD. This suggests that the idea of "thermalization" at an angoricity is able to describe the 
jamming system very well. 
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Thermodynamic analysis of the jamming transition. — So far we have considered 
how the angoricity determines the pressure fluctuations in a jammed packing at a flxed (j). 
The role of the compactivity in the jamming transition can be analyzed in terms of the 
entropy which is easily calculated in the microcanonical ensemble from the density of states. 
Figure m shows S = ln{Q{p, 0)) {Q is the number of states which is the unnormalized version 
of g{r,(j))), which is the non-equilibrium entropy of the system at the given (p, 0) in phase 
space. 

We analyze the non-equilibrium entropy surface S'(ln(0 — 0^), Inp) plotted versus (ln(0 — 
0c), Inp) in Fig. HJand demonstrate that the MD curve (p(0))dyn passes along the maximum 
of the entropy surface constrained by the coupling between p and 0, Eq. Thus the points 
along the entropy surface deflned by (p(0)) correspond to the equilibrium entropy; such a 
curve is superimposed to the entropy surface in Fig. |H Due to the coupling through the 
contact force law, the maximization of entropy is not on p or alone but on a combination 
of both. The entropy S reaches a maximum at the point S'(ln((0)dyn — 0c), In(p)dyn) when we 
move along the direction perpendicular to the jamming curve (p(0))dyn (see the maximization 
direction in Fig. H]). This is a direct veriflcation of the second- law of thermodynamics: the 
dynamical measures maximize the entropy of the system. 

We can use this result to obtain a relation between angoricity and compactivity. We write 
Inp = Inpo + d ln(0 — 0^) where a is the exponent in Eq. (|3]), such that 5'(ln(0 — 0^), Inp) is 
maximized at the MD measures according to the direction of (— sin 6', cos 6') (tan^ = a is the 
slope of the power-law curve in the log — log plot in Fig. H]). Therefore, neither A nor X can 
play the role of the temperature of the system alone, but a combination of both determined 
by entropy maximization satisfying the coupling between stress and strain. Since ^S* = at 
(ln((0)dyn— 0c), In(p)dyn) aloug (— sin6', cos 6*) then (dS/dlnp) cos6'— ((95/(9 ln(0— 0c)) sin6' = 
0. We obtain Cia + ac2/3 = (where Ci = F and C2 = (0 — (f>c)iNVg/(f>'^)) and the relation 
between A and X (see Fig. [13] and Materials and Methods Section Till E^ : 

X = -aA(0-0c)/p0. (8) 
From Eq. (|8]) we obtain that: X oc —(0 — 0c)"^^"~'*'/0 and near 0c: 

X~-(0-0c)^ (9) 
We notice that the compactivity is negative near the jamming transition. A negative tem- 
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perature is a general property of systems with bounded energy like spins [27|: the system 
attains the larger volume (or magnetization in spins) at 0c when X — ?■ 0^ and not X — > +00 
[The bounds 0c < < 1 imply that the jamming point at X — )■ 0~ is "hotter" than 
X — )■ +00. At the same time A — t- 0^ since the pressure vanishes]. 

We conclude that, A and X alone cannot play the role of temperature. Instead, there 
is an actual "jamming temperature" Tj that determines the direction (— sin^,cos^) in the 
log — log plot of Fig. m along the jamming equation of state (see Fig. [T3l) . By maximizing 
the entropy along this direction we obtain Tj as a function of A and X (see Materials and 
Methods Section [Tlll!) : 



= sin4 = ^t4=^^ <PcV-^- (10) 

By the definition of "heat" capacity, we obtain two jamming capacities as the response 
to changes in A and X: 

Cv = dV/dA ~ (0 - 0c)-^ ~ and = dV/dX ~ (0 - 0c)~^ ~ \X\'^I^. (11) 

From Eq. (11 II) . the jamming capacities diverge at the jamming transition as A — i- O"*" 
and X — 7- 0~. However, this result does not imply that the transition is critical since 
from Einstein fluctuation theory applied to pressure and volume j^^l we obtain (we consider 
ks = 1 for simplicity): 

((Ar)2) = A^Cr ~ A^-', and {{AVf) = X^Cy ~ \X\^-'>. (12) 

Thus, the pressure and volume fluctuations near the jamming transition do not diverge, 
but instead vanish as A — )■ 0+ and X — )■ 0~. From a thermodynamical point of view, the 
transition is not of second order due to the lack of critical fluctuations. As a consequence, 
no diverging static correlation length can be found at the jamming point during isotropic 
compression. However, other correlation lengths of dynamic origin may still exist in the 
response of the jammed system to perturbations, such as those imposed by a shear strain 
or in vibrating modes 3, S]- Such a dynamic correlation length would not appear in a 
purely thermodynamic static treatment as developed here. We note though that responses 
to shear can be treated in the present formalism by allowing the inverse angoricity to be 



tensorial 



10| . The intensive jamming temperature Eq. ( ITOj) gives use to a jamming effective 



energy Ej as the extensive variable satisfying Tj = dEj/dS and a full jamming capacity 



Cj ~ (0 — 0c)~^, which also diverges at jamming (see Materials and Methods Section fill ED . 
However, the fluctuations of defined as {{AE^Y) = T|Cj ~ Tj has the same behavior 
as the fluctuations of volume and pressure, vanishing at the jamming transition Tj — )■ O"*" 
[A 0+ in Eq. (^]. 



II. CONCLUSIONS 



We have suggested that the concept of " thermalization " at a compactivity and angoric- 
ity in jammed systems is reasonable by the direct test of ergodicity. The numerical results 
indicate that the full canonical ensemble of pressure and volume describes the observables 
near the jamming transition quite well. From a static thermodynamic viewpoint, the jam- 
ming phase transition does not present critical fluctuations characteristic of second-order 
transitions since the fluctuations of several observables vanish approaching jamming. The 
lack of critical fluctuations is respect to the angoricity and compactivity under isotropic 
compression in the jammed phase — )■ 0^, which does not preclude the existence of critical 
fluctuations when accounting for the full range of fluctuations in the liquid to the jammed 
phase transition from below (pc- Thus, a critical diverging length scale might still appear 
as — )■ c 



29 



30[. Our results suggest an ensemble treatment of the jamming transition. 



One possible analytical route to use this formalism would be to incorporate the coupling 
between volume and coordination number at the particle level found in la| together with 
similar dependence for the stress to solve the partition function at the mean field level. This 
treatment would allow analytical solutions for the observables with the goal of characterizing 
the scaling law near the jamming transition. 



III. MATERIALS AND METHODS 



A. System Information 

The systems used for both, ensemble generation and molecular dynamic simulation, are 
the same. They are composed of 30 spherical particles in a periodic boundary box. The 
particles have same radius r = 5/im and interact via a Hertz normal repulsive force without 
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friction. The interaction is defined as: 

Fn = l ky/\6rY, (13) 

where 6r = (l/2)[2r — \xi — X2I] > is the normal overlap and kn = 4(7/(1 — u) is defined 
in terms of the shear modulus G and the Poisson's ratio u of the material from which the 
grains are made and 6 = 3/2. Here, we use G = 29 GPa and u = 0.2 for spherical particles 
and the density of the particles, p = 2 x 10^ kg/m^. 



B. Ensemble Generation 



In this section, we first explain the method to obtain geometrically distinct minima in 
the PEL. Then we show that the density of the states, ^'(r, 0), does not change significantly 
after sufficient searching time for the configurations. 

For structureless particles possessing no internal orientational and vibrational degrees 
of freedom, at a fixed volume fraction 0, the potential energy is a 3 A^— dimensional function, 
£'(ri, . . . , Tiv), depending on the positions of the N particles. In principle, if all local 
minima corresponding to the mechanically stable configurations of the PEL are obtained, 
the density of states g{r,(j)) can be calculated. Such an exhaustive enumeration of all the 
jammed states requires that not be too large due to computational limits. On the other 
hand, in order to obtain a precise average pressure in the MD simulation, (p)dyru N cannot 
be too small such that boundary effects are minimized. Considering these constraints, we 
choose a 30 particle system. 

In order to enumerate all the jammed states at a given volume fraction 0, we start by 
generating initial unjammed packings (not mechanically stable) performing a Monte Carlo 
(MC) simulation at a high, fixed temperature. The MC part of the method applied to 
the initial packings assumes a flat exploration of the whole PEL. Every MC unjammed 
configuration is in the basin of attraction of a jammed state which is defined as a local 
minimum in the PES with a positive definite Hessian matrix, that is a zero-order saddle. 
In order to find such a minimum, we apply the LBFGS algorithm provided by Nocedal 
and Liu [22]. The procedure is analogous to finding the inherent structures 23| of glassy 



systems. The LBFGS algorithm is also similar to the conjugate gradient method employed 
by O'Hern et al. [4, ll|, but it is computationally more efficient since it does not require 
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the calculation of the Hessian matrix of the system at any time step. The PEL for each 
fixed (f) likely includes millions of geometrically distinct minima by our simulation results. 
Therefore, an exhaustive search of configurations is computationally long. We check that 
the number of found configurations has saturated after sufficient trials such that the density 
of states g{T , 0) has converged to a final shape. 

It is also important to determine if the local minima are distinct. Usually, the eigenvalues 
of the Hessian matrix at each local minimum can be used to distinguish these mechanically 
stable packings. Here, we follow this idea to compare minima to filter the symmetric pack- 
ings. However, instead of calculating the eigenvalues of each packing, which is very time 
consuming, we calculate a function of the distance between any two particles in the pack- 
ing to improve search efficiency (for the LBFGS algorithm, we do not need to calculate the 
Hessian matrix). For each packing, we assign the function Qi for each particle in the system: 

Q^= E *^^'(^)' (14) 

where Vij is the distance between particles i and j, L is the system size and = 30 is the 
number of the particles of the system. We list the Qi for each packing from minimum to 
maximum {Qi}{l < i < N). Since Qi is a higher order nonlinear function, we can assume 
that two packings are the same if they have the same list. The tolerance is defined as: 

^=V ""iV2 ' (^^) 

where Qi and Q[ are the corresponding values from the lists of two packings. 

Figure [5] shows the distributions of the tolerance T for packings at different volume 
fractions. This figure suggests that two packings can be considered the same if T < 10~^, 
which defines the noise level. 

From Fig. [6l we see that after one week of searching, g{r, (p) does not change significantly, 
since the initial packings are generated by a completely random protocol. We also check the 
probability (defined as where A^ncw(^) is the number of new configurations found on 

the ith day and A^totai(0 is the total number of configurations found in i days) of finding new 
mechanically stable states for different searching days. From Fig. U\ we see that, after one 
week searching, the probability of finding new configurations at different volume fractions 
converges, suggesting that enough ensemble packings have been obtained to capture the 
features of g{T, 0). A further test of convergence is obtained below in Fig. [TTl 
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C. MD Generation 



In order to analyze numerical results, we perform MD simulations to obtain Zdyn, Pdyn 
and 0Hvn) which are herein considered real dynamics. The algorithm is described in detail 



in [6|, |16|, |2^. Here, a general description is given: A gas of non-interacting particles at an 
initial volume fraction is generated in a periodically repeated cubic box. Then, an extremely 
slow isotropic compression is applied to the system. The compression rate is Tq = S.Qtg ^, 
where the time is in units of to = -Ra/ p/G. After obtaining a state for which the pressure p 
is a slightly higher than the prefixed pressure we choose, the compression is stopped and the 
system is allowed to relax to mechanical equilibrium following Newton's equations. Then the 
system is compressed and relaxed repeatedly until the system can be mechanically stable at 
the predetermined pressure. To obtain the statistical average of Z^y^ and ^dyn, we repeat the 
simulation to get enough packing samples having statistically independent random initial 
particle positions. Here, 250 independent packings are obtained for each fixed pressure (see 
Fig. E). 



D. Angoricity Calculation 

Since we obtain g{r, 0) and (p)dyn for each volume fraction 0, we can calculate the inverse 
angoricity a by Eq. ([I]). The pressure {p{a,(f)))cns for a given is a function depending on 
a as: 



{p{a, 0))cns = roc ^^ r^-n = • (16) 



/,°°p^(r,0)e-"^dr ^ Epe-"" 

/~^(r,0)e-rdr Ee- 
Figure [9] shows the result of the numerical integration of Eq. ( |T6l) for a particular = 
0.614 as a function of a using the numerically obtained g{r,(j)) from Fig. [H To obtain 
the value of a for this 0, we input the corresponding measure of the pressure obtained 
dynamically (p(0))dyn and obtain the value of a as schematically depicted in Fig. [91 The 
same procedure is followed for every (see Fig. [10]) and the dependence a;(0) is obtained. 
The result is shown in Fig. [JB in the main text. 

We also check the inverse angoricity a(0) using g(r,(j)) for different searching days (see 
Fig. [6]) to ensure the accuracy and convergence to the proper value. From Fig. [Til we can 
see that, after 10 days searching, q;(0) is stable due to the fact that the density of state, 
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g{r, 0), does not change significantly. For volume fraction much larger than 0c, the system's 
input pressure (p(0))dyn reaches the plateau at low a of the function (p(a,0))ens (see Fig. 
[TU]) and the corresponding a(0) becomes much smaller (the angoricity A{(f)) becomes much 
larger), leading to large errors in the value of A as becomes large. This might explain the 
plateau found in A when (0 — 0^) > 2 x 10~^ as shown in Fig. [2J3. 
Using a(0) for each volume fraction, we calculate {Z)ens by: 

/-z,(r,0)e-^dr E^e-r 

the average force (-F)ens by: 

Jo £/(r,0)e "^dr 

where F is the average force for each ensemble packing and the force distribution Pens{F/F) 
by: 

_ /-p(F/F)^(r,0)e-"^dr _ ZPiF/F)e--' 
^"^^ /o°°^7(r,0)e-rdr Ee-"^ ' ^ ^ 

Equations (fT7|) - (fT9|) are then compared with the dynamical measures for a test of ergodicity 
in Fig. |3]in the main text. 



E. Entropy Calculation 

Here we present the calculation of the "jamming temperature" Tj and the corresponding 
jamming "heat" capacity Cj. From the power-law relation p = T/V (x {(p ~ 0c)", we have: 

Inp = Inpo + ctln(0 — 0c), (20) 

where Pq is the constant depending on the system and the slope tan 6' = a. Figure H] indicates 
that the jammed system always remain at the positions of maximal entropy ^S* = in the 
direction (— sin 0, cos 0), perpendicular to the jamming power-law curve. In order to further 
analyze this result, we plot the entropy distribution along the direction (— sin6',cos^) in 
Fig. [121 We see that the entropy of the corresponding jammed states remains at the peak 
of the distributions along (— sin 6', cos 6'), verifying the maximum entropy principle in this 
particular direction. We notice that some deviations are found in the vicinity of 0c. The 
maximization of entropy is not on F or alone, but on a combination of both. This means 
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that the entropy 5'(ln((0)dyn — 0c), ln(j9)dyn) is maximum along the direction of (— sin 6',cos 9) 
and the slope of the entropy along this direction (— sin 6',cos is (see Fig. [13]), that is, 

9S . , dS ^ 

sm9 = — — cos^. (21) 



d\n{(j) — 0c) d\np 
Thus we verify the second law of thermodynamics for jammed systems: 55* = at 

(ln((0)dyn - 0c),ln(p)dyn)- 

By the definition of angoricity A = dV/dS and compactivity X = dV/dS, we have: 

dS dS dS r ci 

Pl^ = ^7^ = ^ = ^^ (22) 



dlnp dp dV A A' 



dS _ ,^dS_ ,^dV I _ C2 



— f23l 

91n(0-0e) ^"'8(1) ^'"d(t)X 02 X X' ^ ' 

where = NVg/V, Ci = T and ca = (0 - 0c)(X1^g/02). By Eq. (ES]) and Eq. (ES]), we can 
simplify Eq. 

The relation between X and A can be obtained then: 

X = —a — A = —a A. (25) 

Ci p(p 

Since we obtain the angoricity A oc {(j)—(j)c)'^ with 7 = 2.5 in the main text and the pressure 
p cc {(f) — 0c)" with a = 1.5 (actually we get 1.65 for the small system size used in the main 
text but the difference can be neglected to simplify). The compactivity X oc —(0 — 0c)^/0. 
We can therefore define the "jamming temperature" Tj as a function of the slope along the 
direction (cos 6*, sin 6*): 

1 Ci . C2 Ci C2. Ci C2 , , 

TTT = -r sm — — cos = cos uia— — — ) = 7: = ~t7 (■^^l 

Tj A X ^ A X^ AsinO XcosO ^ ' 

That is: 

_ y4 sin 6* X cos B sin 6^ , a A , , , . ^ „ , , , . . 

Tj = = = ^A = -j===j, ~ (0 - 0c)"-'^ ~ (0 - 0c). (27) 

Ci C2 i vl + 0. J- 

Furthermore, the "jamming energy" Ej, corresponding to the "jamming temperature" Tj 
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in Eq. (1271) . has the relation as below: 



dE, = T,dS 



, X cos 9. , C2 , , , , , , , A sin ^ Ci , , 

)(-4)dln(0 - 0J + 4dlnp 



C2 X' Ci A 

COS 0d ln((/) — + sin 6'd In p 

(cos 6* + sin 9 tan 6')d ln(0 — cf) 

dln(0 - 0e) 
cos^ 



That is, 



dEj = Va2 + ldln(0 - 0c), 



and 



Ej = (Va2TT) ln(0-0,). 
The jamming capacity Cj can be obtained as: 

dS ^ dS dlnp ^ dS 91n(0 



C,i = Tj— = Tj- — ^ + Tj 



aTj ■'Slnp (9Tj '(91n(0-0c) 5Tj 
Finally, with Eq. (12T]) -( 123|) . the capacity Cj can be calculated: 

_ ci C2 , c^lnp _ 1 + ci glnp 

Since Tj ~ (0 — 0c) and p ~ (0 — 0c)^'^, we obtain Cj ~ (0 — 0c)~"'^- 
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Fig. [T] Ensemble calculations. The density of states g{r,(f)) as a function of internal virial 
r for different volume fraction, 0, ranging from 0.610 to 0.670. The inset is a schematic 
two-dimensional potential energy landscape surface. The jammed states A and B are local 
minima (zero order saddles) in the PES where the external force for each particle is zero 
and the Hessian matrix of the system is positive definite. Our simulation system consists 
of 30 frictionless spherical particles interacting by Hertzian forces with periodic boundary 
conditions. 

Fig. [21 Scaling of pressure and inverse angoricity. (A) The blue Q shows the power-law 
relation for (p)dyn vs (0)dyn — 0c for the 30-particle system. Here, the pressure (p)dyn are 
average values obtained by 250 independent MD simulations. The red Q is the pressure 
used to obtain the inverse angoricity a predicted by Eq. The relatively small system 
size results in large fluctuations of the observables. In order to predict a precise relation 
for the system (A^ = 30), sufficient independent samples of the packings are generated to 
calculate the precise average for observables. We prepare 250 independent packings for each 
to get enough statistical samples to obtain (p)dyn and {Z)dyn by statistical average (see 
Fig. [8]). The inset shows a semi-log plot. (B) The inverse angoricity a as a function of 0-0c. 
We find a power-law relation for system's volume fraction near 0^. The solid line has a 
slope of -2.5. The inset is the angoricity A{= 1/a) vs 0-0c. The plateau observed in A for 
large volume fraction might be related to the finite size of the sample. 

Fig. m Test of ergodicity. (A) The blue Q is the average coordination number {Z)dyn 
obtained by 250 independent MD simulations. The red Q is the coordination number {Z)ens 
calculated by the ensemble for different volume fractions. Agreement between both measures 
supports the concept of ergodicity in the system. (B) The same as (A) but in a log-log plot. 
The blue Q shows the power-law relations for (Z)dyn--^c vs (0)dyn -0c for 30-particle system 
with 0c = 0.6077 and = 5.82. (C) Comparison of (-F)dyn and (-F)ons for different volume 
fractions. (D) The comparison of selected distributions of forces Pdyn{F/F) and Pens{F/F) 
for different volume fractions. 

Fig. m Microcanonical calculations. The entropy surface S'(ln(0 — 0c),lnp). The color 
bar indicates the value of the entropy. The superimposed blue Q is (p(0))dyn from MD 
calculations as in Fig. [2^. The olive arrow line indicates the maximization direction of 
the entropy (— sin^, cos6'). Following this direction, the entropy is maximum at the point 
(ln((0)dyn — 0c), In(p)dyn), Corroborating the maximum entropy principle. 
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FIG. 2: 
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FIG. 5: The distribution of the tolerance T between any two packings at the given (j). From the 
graph, the value of T for which any two different packings are considered to be same is chosen to 
be 10~^, which is above the noise threshold and below the distribution of T. 
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FIG. 6: The distribution of g{T,(p) for 15 days searching (A) at (j) = 0.609, (B) at </> = 0.614, (C) 
at (f) = 0.625. Different color in (A), (B), (C) corresponds to the different day. We find that after 
15 days the distributions have converged. 
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FIG. 7: The probabiUty to find new configurations for different searching day. 



FIG. 8: The cyan Q is ■^dyn and </>dyn for every single packing obtained with MD and the blue Q 
are the average over the single packings for the system which are then shown in the main text of 
the paper. 
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FIG. 9: The numerical integration of Eq. ()16p for = 0.614 is shown as the pink curve. We input 
the (p)dyn (pink Q in the plot) and obtain the corresponding inverse angoricity a. 
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FIG. 11: Calculation of inverse angoricity a for different searching day. 



FIG. 12: The distribution of entropy ^(Inp, ln((/) — (/)c)) along the direction (— sin0,cos0) for 
different jamming ensemble points. The blue Q ^"^^ the entropy for jammed system, which is the 
maximum of S*, verifying the second law of thermodynamics. We notice that some deviations are 
found near 
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FIG. 13: The representation of the maximization analysis SS = along the direction (— sin 9, cos 9) 
for one point in the jamming power-law curve. Here ci = F and C2 = — (f>c){NVg/(tP). 
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